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Abstract 

Assessing solution error continues to be a formidable task when numerically solving practical flow problems. 
Currently, grid refinement is the primary method used for error assessment. The minimum grid spacing requirements 
to achieve design order accuracy for a structured-grid scheme are determined for several simple examples using 
truncation error evaluations on a sequence of meshes. For certain methods and classes of problems, obtaining design 
order may not be sufficient to guarantee low error. Furthermore, some schemes can require much finer meshes to 
obtain design order than would be needed to reduce the error to acceptable levels. Results are then presented from 
realistic problems that further demonstrate the challenges associated with using grid refinement studies to assess 
solution accuracy. 
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1. INTRODUCTION 

A grid convergence study can be used to compute the order of a numerical method. Generally, if the computed 
order matches with the theoretical value, the discrete solution is expected to be a good approximation to the exact 
solution of the equations being solved. This procedure is one of the few methods available to assess the correctness 
of a solution. However, for realistic problems, the cost of obtaining multiple solutions within the asymptotic range 
is often prohibitive. Furthermore, demonstrating the design order for a model problem still does not quantitatively 
reveal the error in the desired solution of a realistic problem. Indeed, some methods can produce highly accurate 
solutions on grids that are well outside of the asymptotic range, and the converse can also be true. Nonetheless, until 
better methods of assessing the error can be developed, grid refinement will remain the primary tool for verification 
and validation of computational solutions. 

Many professional societies 1 ' 2 and journals have adopted guidelines for mesh refinement studies to demonstrate 
that a Computational Fluid Dynamics (CFD) result is a meaningful solution to the governing equations being solved. 
In 1998, Oberkamp and Blottner 3 reviewed the sources of errors in numerical solutions which ranged from boundary 
conditions to inadequate resolution. All of the issues raised in that paper are still important today including the one 
they identified as the largest contributor to solution error - grid resolution. They observed that solutions to the steady 
Navier-Stokes equations often require grids that have an error (relative to the exact solution) as small as 0.1% if the 
expected order property is to be observed in a grid convergence study. They attributed the high accuracy requirement 
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to localized high gradients caused by the nonlinearities in the equations. Most CFD solutions are computed with 
meshes near the limit of what can be afforded, and they rarely produce relative errors as small as 0.1%. Therefore, 
obtaining multiple solutions with such low levels of error required for a grid convergence study is often prohibitive. 

Including time in the refinement of a three-dimensional problem results in a work increase of a factor of r 4 
between any two levels. Here, r is the refinement factor. Typically, three grid levels are used to demonstrate the proper 
relationship between step size and errors. Clearly, the cost would be excessive to perform these calculations unless 
the solution on the coarsest level could be obtained very quickly, which would likely mean that the calculation would 
not be within the domain of convergence of the operators. The simplest procedure for generating a sequence of grids 
is to delete every other point from a finer, structured mesh. This corresponds to a refinement factor of r = 2, but using 
r = V 2 would lead to more manageable overall ratio of effort of 16 for three levels. However, these grids are not 
simple subsets of each other, so they would need to be generated independently, taking care to ensure that they are 
related. 

Another major complication arises for real-world applications involving complex geometries. The formal order of 
accuracy for most methods is only obtained on smooth meshes with smooth solutions. Some important flow features 
involve locally strong gradients that can be very sensitive to the mesh. Block interfaces or overlap regions are required 
for structured grid codes to handle even moderately complex geometries, and these interfaces as well as other boundary 
conditions can introduce additional errors in localized regions. 

For unstructured grid methods, the entire notion of a sequence of refined meshes is difficult to achieve. Nonethe- 
less, researchers have proposed the notion of families 4 of meshes maintaining similar stretching factors that have been 
used with some success in refinement studies. Despite the impediments, mesh refinement is still the primary means of 
assessing the quality of a numerical solution. Before proceeding with some examples of grid convergence studies on 
realistic problems, some simple test cases will be used to examine the resolution requirements to achieve asymptotic 
behavior with several spatial finite-difference operators. 

2. THE ASYMPTOTIC RANGE 

The requirement of being within the asymptotic range to achieve design order accuracy is well known, but what 
kind of resolution that corresponds to is not. Isolating the resolution requirements in a realistic problem is difficult 
because there are usually many competing sources of error. Furthermore, the particular problem solution has a strong 
influence on the results. An analysis in an ideal situation may not directly translate to real-world problems, but the 
estimates should at least provide a lower bound on the resolution requirements. To investigate the spatial accuracy, 
the simplest test - truncation error - involves comparing the exact derivative of a known, analytic function with its 
numerical approximation. The two “trial” functions shown in Figure la are used to test several finite-difference 
operators. The first function is a single-frequency. Cosine function, u\ = 1 /2(cos(7w) + 1), and the second is a 
Gaussian, U 2 = expf-v 2 ). The autospectra in Fig. lb display the energy in each Fourier frequency (normalized by the 
total energy) as a function of the number of wavelengths in the domain. 

Four numerical approximations of the first derivative of the cosine function ui using central-differences are ex- 
amined in Fig. 2: second-, fourth-, and sixth-order operators along with a dispersion-relation-preserving (DRP) 
optimized operator with a 7-point stencil that is formally fourth-order accurate. 5 The plot of the L 2 norm of the error 
versus the number of grid points in Fig. 2a shows that all of the non-optimized operators rapidly approach the design 
order, whereas the optimized scheme transitions to the design order later. When less than 10 points are used to rep- 
resent the Cosine function, the optimized scheme has the lowest error. However, as shown in Fig. 2b, the optimized 
scheme doesn’t approach the design order until there are nearly 30 points within the domain. Indeed, the computed 
order can vary wildly when between 4 and 15 points are used even though the operator is giving a fairly accurate 
representation of the derivative. In general, the higher the order of the operator, the more points are needed to achieve 
the design order, although all of the standard operators are very close to design order with 15 points. 

A similar analysis for the Gaussian function is presented in Fig 3 where the results from a uniform grid are 
compared with those from a stretched grid. The clustered mesh points were generated using x = -4 + 4(1+ sinh(3(/ ; - 
1/2))/ sinh(3/2)) where 0 < £ < 1 in the computational space. The derivatives are then d/dx = (1 /X()d/d^. 

The observed error in the derivative of 112 in Fig- 3a is oscillatory as the grid count is increased up to around 12, 
then begins to decrease with mesh size. Clearly, the resolution must be sufficient to resolve more than the first Fourier 
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(b) Autospectra of trial functions. 


Figure 1: Trial functions used in grid convergence studies. 



Figure 2: Grid convergence study for the derivative of Cosine function u\. 
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frequency to get a reasonable result. Clustering the grid around the region of high gradients does reduce the error, 
but only once the waveform is sufficiently resolved. The computed orders in Fig. 3b are also highly oscillatory until 
there are around 15 points in the domain. The order of the optimized scheme oscillates until there are approximately 
50 points representing the Gaussian. To achieve design order with the optimized scheme would require over 100 
points to represent the Gaussian. Even though most of the energy is in the first frequency, the higher frequencies must 
be resolved to achieve design order. The second-order method approaches design order the fastest, but it is also the 
least accurate. By examining Figs. 3a and 3b together, the lack of direct correlation between accuracy and achieving 
design order is apparent. In most applications, an error between 10 2 and 1 0 1 would be acceptable. The higher-order 
methods achieve this well before approaching the asymptotic regime. However, the second-order method requires 30 
to 100 points for this level of accuracy which means that it is operating within the asymptotic regime. Certainly in 
this example, the notion of using grid convergence to assess accuracy is only appropriate for low-order methods. 



Figure 3: Grid convergence study for the derivative of Gaussian function U 2 using a uniform grid (dashed lines) and 
sinh stretching (solid lines). 

The solution of most problems of interest will involve the approximation of higher derivatives. Fig. 4 shows the 
results for the second derivative of U 2 using two applications of the first-derivative operators as well as by second- 
derivative operators. The second-derivative operators are more accurate, but both methods achieve design order at 
about the same rate. Actually, the behavior of the computed order for the second derivative is not much different than 
that obtained for the first derivative. 

Although this exercise does give an indication of the resolution requirements to achieve design order, real world 
problems are much more complex. The Navier-Stokes equations involve combinations of multi-dimensional first- 
and second-derivative spatial operators dealing with drastically varying scales throughout the domain. Furthermore, 
the nonlinearities in the equations can result in changing scales as the mesh is refined. Boundary closure, mesh 
smoothness, and many other factors will play into the actual order that is observed in a solution. The points-per- 
wavelength estimates seen in these simple examples were quite high for even a smooth Gaussian distribution, so 
observing design order in realistic problems is likely to be challenging. 

3. GRID CONVERGENCE EXAMPLES 

According to the ASME guide on Verification and Validation, 1 “systematic grid refinement is the cornerstone 
of the verification process for either codes and solutions.” Code verification confirms that a model has been coded 
correctly. The method of manufactured solutions 6 ' 7 is a procedure for generating a forcing function to the governing 
equations such that an arbitrary analytically defined manufactured solution satisfies the governing equations including 
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Figure 4: Grid convergence study for the second derivative of the Gaussian function ih using sinh stretching. Solid 
lines for second derivative operators and dashed lines for two applications of a first derivative. 


the forcing function. As the grid is refined, the numerical solution can be compared against the manufactured solution 
to ascertain whether the errors follow the expected asymptotic behavior. Hence, the errors can be evaluated directly. 
Multiple manufactured solutions are typically needed to exercise every term in the governing equations. The costs of 
the tests can be reduced dramatically by examining the truncation error on subsets of the entire grid. A truncation error 
test only requires a single evaluation of the numerical residual to compare with the exact residual of the manufactured 
solution. Furthermore, the problem can be downscaled by choosing focal points within the mesh and only evaluate 
the residual within windows around the focal points that shrink as the mesh is refined. Thomas et al. s have shown 
that the estimate of the order from truncation error analyses on downscaled domains are good indicators of the order 
that would be obtained using the full grid. 

In solution verification, the error is estimated using Richardson extrapolation. 9 In general, three grid levels are 
needed to ascertain the order of the scheme and estimate the error. The ASME recommends a refinement factor 
between grid levels of at least 1.3. However, Richardson extrapolation assumes that the solution on all three levels 
is within the asymptotic range of the convergence of the operators. For realistic problems, this level of refinement is 
difficult to achieve. Vassberg and Jameson 10 found inconsistent and somewhat disappointing convergence properties 
of three CFD codes for steady, inviscid flow over a NACA0012 airfoil. Two examples are included in this section to 
demonstrate typical results from grid refinement in unsteady, viscous problems. 

3.1. Tandem Cylinders 

A tandem cylinder arrangement with a centroid to centroid spacing of 3.7 diameters was investigated numerically 
in reference 1 1 . The inline cylinders are a prototype for the type of bluff-body interaction that occurs around landing 
gear. The Mach number was 0.166, and the Reynolds number based on the diameter was 166 thousand. However, the 
simulations were intended to match with experimental data collected with transition strips installed on the cylinders. 
The experimental surface pressure distributions resemble those obtained at Reynolds numbers above 4 million. 

Figure 6 shows a planar cut of the coarsest mesh used in this study. The physics in the problem is fairly compli- 
cated, and the grid needs to be sufficiently fine to resolve all of the important scales demonstrated in Fig. 5. Although 
the geometry is relatively simple, the block-structured grid still includes block interfaces that are not one-to-one to 
reduce the overall grid count. Furthermore, the grid is fairly distorted in the region where the topology switches from 
circular to rectangular around each cylinder. The calculations were performed using the CFL3D 12 Navier-Stokes flow 
solver, which is a second-order method in both space and time. A k-aj SST based hybrid RANS turbulence model 
was employed with the production term turned off at a fixed distance from the cylinders roughly corresponding to the 
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thickness of the boundary layers. 11 Hence, the method basically becomes an implicit Large-Eddy-Simulation (LES) 
method outside of a fixed region enclosing the boundary layers. Although traditional LES methods have a varying 
filter width as the mesh is refined, implicit LES relies upon the artificial dissipation in the scheme to provide the 
dissipation of high wave numbers. Hence, one can still expect to observe an order property from scheme as the mesh 
is refined. Nonetheless, the mesh must be sufficiently fine to resolve all of the dominant scales in the flow. The case 
was run fully turbulent to best simulate the experimental configuration with transition strips. 



Figure 5: Density contours around the tandem cylinders from fine grid solution. 


>> 0 


Calculations were performed on grid levels with 
total point counts of 1.35, 10.1, and 78.3 million. 

The succession of meshes were generated by elim- 
inating every other point in all directions from the 
mesh one level finer. The time step was varied on the 
finest mesh to assess the required time step, then the 
grid refinement calculations were performed with the 
same time step to isolate the effects of the mesh. The 
variation of the lift and drag along with their pertur- 
bations (rms of the perturbations are denoted by ') 
with mesh size is presented in Figure 7. The total 
number of points is given by N, so N~ 1/3 gives an 
indication of the spacing in each of the coordinate 
directions. As CFL3D is a second-order method, 
the variation of lift and drag within the asymptotic 
regime should be linear when plotted against N~ 2/3 . 

Clearly, only the drag on the downstream is even 
close to linear on this scale, and the randomness in 
the trends is a good indication that this result is fortu- 
itous. Indeed, the computed order properties in table 

1 are not near the expected value of 2. Furthermore, the order cannot be computed for the lift because the trend is not 
continuous between the three meshes. The time-averge of the lift should be zero, so the calculated mean lift may give 
an indication of the error caused by the limited temporal sampling length rather than just inaccuracies caused by grid 
resolution. 

The grid convergence index (GCI) 6 is also presented in table 1. A conservative value of the factor of safety 
F s — 3.0 is used in the formula 



Figure 6: Planar view of cylinder grid. 


GCI fie = F * 


Kr p 21 - 1) 


where 0 is the solution variable of interest and the subscript indicates the grid level, with 1 being the finest level. The 
ratio of the grid spacings on successive levels is denoted by r, and p is the observed order of the scheme. The GCI 
(expressed as a % of the solution value) is a semi-empirical estimate of the 95% confidence interval limits. The GCI is 
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fairly low for the downstream cylinder where the observed order is relatively large. However, the lift and drag on the 
downstream cylinder is merely a reflection of what is shed from the upstream cylinder. If the solution on the upstream 
cylinder is not well converged, it cannot be on the downstream cylinder. Hence, we are observing anomalous behavior 
in solutions that do not demonstrate asymptotic convergence. 



Figure 7: Results of grid convergence study on the tandem cylinder configuration. 


Table 1: Computed order property and GCI for the tandem cylinder configuration. 



C L 

c D 

C L 

c 

Cylinder 

order 

GCI 

order 

GCI 

order 

GCI 

order 

GCI 

Upstream 

- 

- 

0.175 

22.2 

0.115 

1070 

0.187 

533 

Downstream 

- 

- 

1.85 

2.66 

2.44 

3.39 

3.72 

0.417 


3.2. Nose Landing Gear 

A family of three unstructured meshes has been generated around a Gulfstream G550 nose landing gear model. 13 - 14 
The FUN3D, 15,16 second-order accurate, unstructured grid solver was used to perform the simulations. The calcu- 
lations were performed with the HRLES 17 hybrid RANS turbulence model. Away from walls, the HRLES model 
transitions to a traditional LES model with a filter width that is dependent on the grid spacing. Hence, additional 
energy at higher frequencies is included in the simulation as the grid is refined which will limit the ability to see grid 
convergence until all physically important scales are resolved. The Mach number was 0. 166 and the Reynolds number 
based on the main strut was 73 thousand. The calculations were run with the turbulence model in a fully turbulent 
mode. Completely tetrahedral meshes were created along with a set with hexagons in the near-wall regions. Previous 
studies have indicated that mixed-element meshes with hexagons near curved walls may be more accurate than the 
fully tetrahedral meshes. 18,19 The sequence of meshes has a refinement factor of V2. The stretching employed in the 
advancing layers is similar for all meshes in a family, but there is some ambiguity in viscous regions where multiple 
advancing fronts are employed. A time step study was performed on the finest mesh, and the grid refinement was 
performed with the same time step to isolate the effects of the mesh. 

Even though the model geometry is not quite as complex as the full-scale flight article, the generation of the 
unstructured grids was difficult. A planar crinkle-cut through the mesh is shown in Fig. 8a. On the coarsest mesh 
of 9 million nodes, the contours of the instantaneous Q criterion in Figure 8b show that there is considerable vortical 
structure in the wake of the gear. A careful examination of the solution on even the finest mesh of 71 million nodes 
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shows that some thin shear layers are not sufficiently resolved. Hence, one would not expect to see the design order 
of accuracy in the grid convergence study. Nonetheless, the convergence appears reasonable for the full tetrahedral 
meshes as shown in Figure 9. The streamwise direction is x, and z is along the axis of the main strut. The convergence 
on the mixed-element meshes is inconsistent, so the order cannot be calculated. Table 2 reveals that the computed 
order for the fully tetrahedral meshes is close to 2. However, this result is almost undoubtedly fortuitous. The cost of 
obtaining a solution on the finest mesh is expensive, and, although the uniform refinement does improve the solution, 
many more levels would be needed to fully resolve some of the flow features. Although one may resort to manually 
specifying local grid refinements, the grid generation task becomes a tedious endeavor, and many of the advantages 
of resorting to an unstructured grid are lost. 


Table 2: Computed order property and GCI for the nose landing gear. 



c x 

Q 

C' 

C 'z 

Grid 

order 

GCI 

order 

GCI 

order 

GCI 

order 

GCI 

Tetrahedra 

2.42 

0.789 

2.19 

5.00 

2.37 

61.3 

3.97 

35.9 

Mixed 

- 

- 

- 

- 

6.89 

2.59 

9.08 

0.919 




(b) Instantaneous Q criterion. 


(a) Planar view of the grid. 


Figure 8: Nose landing-gear grid and results colored by perturbation pressure from the 9 million node mesh. 


4. CONCLUDING REMARKS 

For the structured grid in the tandem cylinder example, the fine grid was generated to have approximately 50 points 
within the boundary layer, so there are only around 12 points to represent the boundary layer variations on the coarsest 
mesh. When using a factor of 2 refinement between levels, the finest mesh would have to be extremely fine to keep all 
of the coarser levels in the asymptotic range. Non-integer refinement must be performed within the grid generation 
process, and grid generation software is only beginning to include this feature. However, even in the example using 
non-integer grid refinement in the context of unstructured grids, the solutions were outside of the asymptotic regime. 
Furthermore, it is questionable whether even the finest meshes in use today produce solutions within the asymptotic 
regime. Hence, solution verification through grid refinement is problematic. Indeed, the conclusion of the 3rd Lisbon 
Workshop on Code Verification, Solution Verification and Validation 20 is that it is safest to assume that data are 
not within the asymptotic range. Nonetheless, the assessment of the accuracy of a numerical solution is extremely 
important. In situations where the error cannot be assessed via grid refinement, alternative methods must be developed 
and used. Currently, comparisons against benchmark experiments and collective comparisons between many different 
codes are techniques that are often employed to develop confidence in predictive capabilities. However, the confidence 
only extends to problems (with grids and physics) that are very similar to those already investigated. Furthermore, no 
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Figure 9: Results of grid convergence study on the nose landing-gear configuration, 
real error estimate is available from such approaches. Such studies are certainly valuable, but more general techniques 
that can be applied to realistic problems are sorely needed. Roy 21 has suggested that mesh adaptation may enable not 
only better solutions but better error estimates because solutions will be driven into the asymptotic regime. However, 
general adaptation algorithms still need development for steady problems. Adjoint methods can be used to guide 
adaption and directly give an error estimate for quantities of interests, 22 but they are difficult to implement and have 
not been demonstrated to be effective within boundary layers. An added complication for aeroacoustics is the unsteady 
nature of the problem. To be effective, the grid would need to adapt in time to the solution which would require that 
the adaption process be efficient and sensitive to the changing nature of the solution. 

In the context of code verification where the cost of the simulations can be controlled and the exact solution is 
known analytically, using the computed order to ascertain the correctness of coding is sensible. However, one rarely 
achieves design order in realistic problems, and the current rate of increase in computer power alone is unlikely to 
enable it in the near future. Grid convergence studies are still useful to demonstrate the sensitivity of the solution to the 
mesh, but they do not always give a complete indication of accuracy. The challenge is to perform the best simulations 
possible with the resources available while performing checks to minimize the likelihood of being misguided by 
erroneous results. Tests with minimal expense are most likely to be performed. Applying some code verification 
methodology, such as a truncation error analysis of the residual on subsets of the realistic mesh, are feasible today. 
However, generating the sequence of grids and finding a manufactured solution that resembles the actual solution is 
still difficult. Furthermore, a convergence study of the truncation error on irregular meshes is not always representative 
of the discretization error. 23 A more desirable option would be to eliminate the need for a sequence of meshes by 
ascertaining something meaningful from the absolute level of the truncation error on a single mesh. Bringing such 
ideas to fruition are needed to meet the challenge of finding practical yet effective ways to quantify and control the 
errors in unsteady simulations around complex geometries. 
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